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The theory of phase stability in the Ni-Au alloy system is a popular topic 
due to the large size mismatch between Ni and Au, which makes the effects of 
atomic relaxation critical, and also the fact that Ni-Au exhibits a phase sepa- 
ration tendency at low temperatures, but measurements at high-temperature 

show an ordering-type short-range order. We have clarified the wide dispar- 
ts . 

ity which exists in the previously calculated values of mixing energies and 
thermodynamic properties by computing "state-of-the-art" energetics (full- 

cn 

' potential, fully-relaxed LDA total energies) combined with "state-of-the-art" 

(N ' 

statistics (k-space cluster expansion with Monte Carlo simulations) for the Ni- 

I Au system. We find: (i) LDA provides accurate mixing energies of disordered 

-t— > ' 

Nii-zAuc alloys (A_£f m i x ^5 +100 meV/atom) provided that both atomic re- 
laxation (a ~100 meV/atom effect) and short-range order (~25 meV/atom) 
q I are taken into account properly, (ii) Previous studies using empirical po 

tentials or approximated LDA methods often underestimate the formation 
^ ' energy of ordered compounds, and hence also underestimate the mixing en- 

ergy of random alloys, (iii) Measured values of the total entropy of mixing 
combined with calculated values of the configurational entropy demonstrate 
that the non- configurational entropy in Ni-Au is large, and leads to a signif- 
icant reduction in miscibillity gap temperature, (iv) The calculated short- 
range order agrees well with measurements, and both predict ordering in the 



disordered phase, (v) Consequently, using inverse Monte Carlo to extract in- 
teraction energies from the measured/calculated short-range order in Ni-Au 
would result in interactions which would produce ordering-type mixing ener- 
gies, in contradiction with both experimental measurements and precise LDA 
calculations. 
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I. INTRODUCTION 



The Ni-Au alloy system is physically interesting because, on one hand it exhibits a phase 
separation tendency at low temperatures and positive mixing enthalpies |l| and, on the 
other hand, an ordering-type short-range order (SRO) at high temperatures. Also, the 
fee Ni and Au constituents possess a large lattice-mismatch (Aa/a ~ 15%), thus making 
this system a critical test for any alloy phase stability theory hoping to capture the effects of 
atomic relaxation. Important early experimental and theoretical work on this alloy includes 
the work of Moss et al. OM, Cohen et al. |5],|2||]], and and Cook and de Fontaine M. The 
coexistence of phase separation (at low T) with short-range ordering (at high T) in the 
same alloy system might have been naively construed to imply a change from repulsive 
("ferromagnetic") interactions at low T to attractive ( "anti- ferromagnetic" ) interactions at 
higher T. The change would have been surprising, given that no electronic, magnetic, or 
structural change is observed in this temperature range. The answer to this puzzle was 
given by Lu and Zunger: || The excess energy for a disordered Mi^Au^ alloy or an ordered 
compound of type a is given by: 

AH a = K« q ) - [(1 - x)E m (a™) + xE Au (a£)}, (1) 

and may be written [[J AH = Ae + AEy-p, where Ae is the constant-volume, "spin-flip" 
energy required to create o out of Ni and Au, each already prepared at the alloy lattice 
constant a a , and AEw is the volume deformation energy required to hydrostatically de- 
form Ni and Au from their respective a^l 1 and a^ u to a*. In Ref. |§, it is demonstrated 
that SRO is determined by the constant volume energy change Ae, which is negative (order- 
ing, or "anti-ferromagnetic") in Ni-Au, indicating an ordering tendency of SRO. However, 
A_Evd = G(x) is large and positive, making AH > 0. And, since long-range order is de- 
termined by AH, Ni-Au shows phase-separating ("ferromagnetic") long-range order. This 
analysis leads to two unexpected conclusions: First, that the time-honored Ising-like repre- 
sentation of alloy thermodynamics which includes only "spin-flip" energies of the Ae type, 
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but ignores the elastic energy G(x) will fail in explaining basic stability trends for systems 
such as Ni-Au. Second, since measurements or calculations of the SRO are insensitive to 
physical effects (i.e., elastic deformation A£Vd) that control measurements/calculations of 
mixing enthalpies AH, the often-used practice [I(J of "inverting" the SRO profile to extract 
interaction energies that are then used to predict mixing enthalpies is fundamentally flawed. 
Specifically, inversion of the SRO of Ni-Au will produce ordering-like interaction energies 
which, when used to calculate mixing enthalpies will produce (ordering-like) negative values, 
while the measured ones are strongly positive. 0,0] 

For these and other reasons, the theory of phase stability in Ni-Au has recently be- 
come quite popular |8|,p^, p!6| , p^ - [T5| , p!9| , p!8| , p!7| (Table |). These calculations are distinguished 



by the methods used for (i) energetics (T=0 K) and (ii) statistics (T 7^ 0). Energy cal- 
culations (T=0 K) for this system have been performed by a wide variety of techniques: 
First-principles calculations, both full-potential (FLAPW || and FLMTO |12|) and atomic- 



sphere-approximation (LMTO [p~3^p~5|] and ASW |16[]), as well as semi-empirical (EAM [ jP7| ) 
and empirical potentials ]IB|JT^ , p!P| . There are significant variations in the computed ener- 
getics (Table |). Statistics have been applied to these calculations using cluster expansions 
(CE) such as e-G ||, Connolly- Williams Pm , and second-order expansions. pi| 
The purpose of this paper is thus three-fold: 

First, we would like to clarify the conflicting energetic and statistical results (Table |) by 
computing "state-of-the-art" energetics for Ni-Au alloys (full-potential LAPW total energies 
including full atomic relaxation) combined with "state-of-the-art" statistics (a k-space CE 
|33| with Monte Carlo simulations). These computations will clarify whether the better 
agreement with experimental AH obtained by approximated methods (e.g., empirical and 
semi-empirical potentials, as well as atomic-sphere-approximation methods) relative to full 
LDA methods is fundamental or accidental. 

Second, we would like to address the issue of why the calculated miscibility gap tempera- 
tures are often much too high compared with the experimentally assessed phase diagram [|TJ . 
In Table |, one can see a fixed ratio between calculated miscibility gap temperatures Tmg 



and the calculated AH m \ x . In fact, all previous calculations (except the EAM calculations of 
Asta and Foiles [|T7|]) very nearly follow the ratio obtained using mean-field configurational 
entropy: &bTmg/ AH m i* = 2. However, the experimental value of this ratio is 1.2. We will 
examine this apparent discrepancy between experimental Aif m i X and Tug below. 

Third, we would like to examine the SRO in Ni-Au and discuss the implications of this 
SRO on "inverse" techniques, mentioned above, for calculating phase stability in alloys. We 
will offer a challenge to practitioners of the inverse Monte Carlo method. 

II. CHECKING ORDERED COMPOUND FORMATION ENERGIES 

Table | summarizes the previous calculations on the mixing enthalpies of random Ni-Au 
alloys. The wide discrepancy between calculated values of AH mix (48-170 meV/atom) is 
apparent from this table. Since mixing enthalpies Aif mix of random alloys can be expressed 
[see, e.g., Eq. (3b) in Ref. ||22|| 1 as a linear combination of formation enthalpies AHf(o~) of 
certain ordered compounds {cr}, the discrepancies in AH mix must reflect discrepancies in 
AHf{cr). But formation enthalpies of small-unit-cell ordered compounds can be computed 
accurately and reliably via full-potential fully-relaxed LDA methods. Our strategy will thus 
be to trace the source of the discrepancy in AH mix to the values of formation energies of 
various Ni p Au g ordered compounds, as shown in Table |TJ. Examining this table leads to 
several interesting points regarding the energetics in Ni-Au: 

A. FLMTO vs. ASA methods (LMTO, ASW) 

In comparing the full-potential LMTO Jl2] to LMTO-ASA |15] calculations, one can 
see significant and strongly configuration-dependent discrepancies, even when considering 
unrelaxed configurations. For example, the Z2 structure [a Ni 2 Au2 (001) superlattice] has 
an unrelaxed formation energy which is nearly 100 meV/atom lower in the LMTO-ASA 
calculation than in the full-potential LMTO one. Thus, the ASA-based calculations (LMTO, 
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ASW) in the Ni-Au system cannot be trusted for the kind of quantitative energetics required 
in phase stability studies. 



B. Harmonic vs. anharmonic relaxation 

In a large lattice-mismatched system like Ni-Au, the effects of atomic relaxation are 
likely to be crucial. Although straightforward, fully relaxing all of the cell-internal and 
cell-external degrees of freedom can be computer intensive. One alternative to full atomic 
relaxation (using quantum mechanical forces and total energy minimization) which has been 
used in Ni-Au || is to use continuum elasticity theory |24] to find the relaxed geometry, with 



a subsequent LDA calculation with this geometry to find the relaxed energetics. Continuum 
elasticity theory can be used as a relaxation model by realizing that many ordering Ni p Au 9 
compounds can be described as "superlattices" along some special orientations k. Contin- 
uum elasticity then provides the equilibrium interlayer spacing c eq along k as a function of 
the externally-fixed perpendicular lattice constant a± as the minimum of the epitaxial strain 
energy due to the external constraint: 

c cq (k,a ± ) = a'J + [2 - 3q W (a ± ,k)[a^ - a J + ... 

(2) 

. AE%(a ± ,k) 
q{a±,k) = (3) 

where and are the equilibrium energy and lattice constant of the cubic material 
A. A-Egpj is the energy of the alloy constituent subject to the biaxial constraint that the 
lattice constant perpendicular to k is externally fixed to be a±. Ai?buik(a±) is simply the 
deformation energy change upon hydrostatically distorting the material from a eq to a±. The 
central quantity in these elastic calculations is the "strain reduction factor" q(a±,k). In 
continuum elasticity theories, q(a±, k) is given by 

q(a ± , k) = 1 - B/[C n + j(a ± , k)A]. (4) 

where 



A = C44 - (Cu - C 12 )/2 (5) 

is the elastic anisotropy, B = (Cu + 2Ci2)/3 is the bulk modulus, and C^- are elastic 
constants. In the harmonic approximation, q(a±, k) is further assumed to be a^-independent, 
and 7 harm '(fc) is the following geometric function for the direction k = (l,m,n): 

4(l 2 m 2 + m 2 n 2 + n 2 l 2 ) 



7harm.(^ m > n ) 



(l 2 + m 2 + n 2 ) 2 
4 2 

-V47T [K (l, m, n) i= k a{1, rn, n)} (6) 

v 21 



where Kl are the Kubic harmonics of angular momentum L. 

Using Eqs. (@)-(@) thus provides predicted relaxed geometries c eq (k,a±) for alloy com- 
pounds (e.g., the Z2 structure) given the elastic constants and . Indeed, these equations 
have been routinely used (see review in Ref. [^5|) to predict the distortion c cq — a cq of films 
grown epitaxially on a substrate with lateral lattice constant aj_. Comparison to LDA cal- 
culations |26| shows that for semiconductors with lattice mismatch (a eq — a±)/a± ^ 7%, the 
harmonic expressions (§)-(§) work very well down to a monolayer thickness. However, we 
find that for noble- and transition-metal alloys with a much larger lattice mismatch (e.g., 
Ni-Au, Cu-Au with Aa/a = 15%, 12%, respectively), anharmonic corrections are impor- 
tant. As we will see below in Sec. |III C| , this is manifested by the fact that 7(a_i_, k) of Eq. 



(|) now has additional terms to those appearing in the harmonic form of Eq. @. These 
anharmonic terms in 'j(k) lead via Eq. (H) to corrections to q(a±, k), and consequently via 
Eq. @ to the relaxation of the lattice constant c eq (k). Indeed, using the same FLAPW as 
Ref. ||, but minimizing the total energy quantum- mechanically ("Fully relaxed" in Table 
|1D rather than via the harmonic expression of Eq. (|J) ("Partially relaxed" in Table ||), 
we find a lower-energy relaxation for Z2: The LDA energy minimization gives AH(Z2) = 
+70.2 meV/atom, while LDA with harmonic relaxation gives +124.3. For other structures, 
the effect is much lower. Nevertheless, anharmonic relaxation in Ni-Au alloys is large and 
cannot be neglected. 
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C. Empirical methods: Getting the right AH m i x (x,T) for the wrong reason? 

We see from Table | that the methods that use empirical evaluations of Aif mix (l/2, oo) 



15| , p^Jlq , p!2|JT7| produce results that are lower, and thus closer to the measured 
AH mix (l/2, 1150) than methods that use converged, full potential, fully relaxed approaches 
(i.e., the present work and Refs. P,|T2"||). Since there is a proportionality between AH m[x 
and AHf(cx), we thus surmise that the empirical methods will produce formation energies 
AHf(a) of ordered compounds that are lower than the LDA results for such systems. In- 
deed, Table [TJ shows the formation energies of two of the empirical potential methods. By 
comparing these numbers to full-potential LDA energies, one can see that the empirical 
potentials systematically underestimate the formation energies of ordered compounds. Since 
the LDA method is expected to reproduce formation enthalpies of small-unit-cell ordered 
structures rather accurately, and since FLAPW gives a precise representation of the LDA, 
we think that the underestimation of FLAPW energies by the empirical methods is a rather 
serious limitation of these methods. The EAM of Ref. [|17| was fit to the unrelaxed FLAPW 
calculations of Ref. ||, and thus reproduces these energies fairly well (except for the Z2 
structure). However, the EAM severely overestimates the energetic effect of relaxation, and 
hence produces relaxed formation energies which are much lower than LDA, and in some 
cases are even negative. |37| It would be desirable to see more formation energies of ordered 
compounds from the empirical methods to determine test the expectation of underestimation 
of AH f (cr) relative to LDA. 

In summary, the reason that empirical methods agree with measured random-alloy mix- 
ing enthalpy better than LDA methods do is a systematic underestimation by the empirical 
methods of even the ordered compound energies. 

III. PRESENT CALCULATIONS - FLAPW WITH K-SPACE CLUSTER 

EXPANSION 



S 



A. FLAPW calculations of ordered compounds 



We have performed first-principles full-potential LAPW |2q] calculations for pure Ni, 
pure Au, and a large number (31) of fcc-based Ni-Au compounds in order to construct 
an accurate cluster expansion. The total energy of each compound is fully minimized with 
respect to volume, cell-internal, and cell-external |^SJ coordinates using quantum-mechanical 
forces. We have used the exchange correlation of Wigner |RU| . The muffin-tin radii are chosen 



to be 2.2 a.u for Ni and 2.4 a.u. for Au. Brillouin-zone integrations are performed using the 
equivalent k-point sampling method, |JT with the k-points for each compound all mapping 



into the same 60 special k-points for the fee structure. This mapping guarantees that the 
total energy per atom of an elemental metal calculated either with the fee unit cell or with 
a lower symmetry (e.g., any of the compounds) are identical. All calculations performed are 
non-magnetic. (The spin polarization energy difference between ferro- and non-magnetic fee 
Ni was calculated and found to be -50 meV/atom.) 



The 31 calculated LAPW formation energies are given in Table |T|. Both relaxed and 
unrelaxed (total energy minimized with respect to volume, but with cell-internal and cell- 
external coordinates held fixed at ideal fee positions) formation energies are shown. The 
nomenclature of the compounds studied is the same as given in |^2] . Many of the compounds 
considered can be described as Ni p Au ? "superlattices" along a particular orientation k: 



NiiAui : 


[100], 


[Ill], 




Ni 2 Aui : 


[100], 


[011], 


[111], 


NiiAu 2 : 


[100], 


[011], 


[111], 


Ni 3 Aui : 


[100], 


[011], 


[201], [111], [311], 


NiiAu 3 : 


[100], 


[011], 


[201], [111], [311], 


Ni 2 Au 2 : 


[100], 


[011], 


[201], [111], [311], 


Ni 3 Au 3 : 


[100], 






Ni 2 Au 3 : 


[100]. 







(7) 



We also calculated the energies of six other structures: Ll 2 (Ni 3 Au and NIA.U3), -D7(Ni 7 Au 
and NiAu 7 ), and two 8-atom "special quasi-random structures" P2[ , SQS14 a (Ni 6 Au 2 and 
Ni 2 Au6). In addition the Ni p Au g long-period superlattice limits (p,q — > 00) needed in the 
construction of the k-space cluster expansion (see below) were computed for six principle 
directions: [100], [Oil], [201], [111], [311], and [221]. The numerical error of the LAPW 
calculations of AHf is estimated to be ~10 meV/atom or less. 

B. k-space cluster expansion 

The Ni-Au formation energies AH a for structures a are then mapped onto a cluster 
expansion using the k-space formulation of Laks et al. |$3f Rather than a cluster expansion 
of AH a , we will expand with respect to a reference energy: 

AE CE (cr) = AH LDA (a) - E rei (8) 

We will separate the CE into two parts: (i) the terms corresponding to pair interactions with 
arbitrary separation will be conveniently summed using the reciprocal-space concentration- 
wave formalism, and (ii) all terms but the pairs will be cast in real-space: 

AE CE (a) = £ J(k)\S(k, a)| 2 + ]T D f J f U f (a). (9) 

k / 

The first summation includes all pair figures and the second summation includes only non- 
pair figures. In the reciprocal-space summation in Eq. @, J(k) and S(k,a) are the lattice 
Fourier transforms of the real-space pair interactions and spin-occupation variables, and 
Si, respectively, and the spin-occupation variables take the value S% = — 1(+1) is the atom 
at site % is Ni(Au). The function J(k) is required to be a smooth function by minimizing 
the integral of the gradient of J(k). The real-space summation of Eq. (|9|) is over /, the 
symmetry- distinct non-pair figures (points, triplets, etc.), Df is the number of figures per 
lattice site, Jf is the Ising-like interaction for the figure /, and 11/ is a product of the 
variables Si over all sites of the figure /, averaged over all symmetry equivalent figures of 
lattice sites. 
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The reference energy of Eq. is chosen to contain infinite-range real-space elastic 
interaction terms. Subtracting these long-range terms from AH^f )A before cluster expanding 
removes the k — > singularity, and thus significantly enhances the convergence of the CE. 
p3| The form used for E ie f is 

£ ref = Ax( l_ x) E AE cs(kx)\S(k,a)\ 2 (10) 

where AE^? s (k,x) is the equilibrium constituent strain energy, defined as the energy change 
when the bulk solids Ni and Au are deformed from their equilibrium cubic lattice constants 
a jvj and olau to a common lattice constant a±_ in the direction perpendicular to k. AE^? s (k, x) 
can thus be written as the minimum of the following expression with respect to a_i_: 

AE e c l(k,x) = (l-x)q N \a ± rk)AEZM + xq Au (a L rk)AE^ k (a ± ). (11) 

where q^ x \a±, k) is given by Eq. (fj). 

The final expression used for the formation energy of any configuration a is then 

AH(a) = £ J(k)|S(k, ^)| 2 + E Df Jf n/M 

k / 

+^hx) ? Ai ^' x) |5(k ' a) |2 (12) 

The following inpui is needed to construct this Hamiltonian for Ni-Au: (i) the formation 
energies of a set of ordered fcc-based compounds (required to fit the values of J(k) and Jf), 
and (ii) the epitaxial energies of fee Ni and fee Au (required to compute the anharmonic 
values of AE^(k,x)). The output is a Hamiltonian [Eq. (0)] which (i) predicts the energy 
of any fcc-based configuration (i.e., not only ordered compounds) even 1000-atom cells or 
larger, (ii) possesses the accuracy of fully-relaxed, full-potential LDA energetics, and (iii) is 
sufficiently simple to evaluate that it can be used in Monte Carlo simulations, and thereby 
extend LDA accuracy to finite temperatures. 
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C. Anharmonic calculation of constituent strain 



Laks et al. |33|] demonstrated that the calculation of AE^ s (k,x) of Eq. flTT|) is signifi- 
cantly simplified if one uses harmonic continuum elasticity theory [i.e., insert Eqs. 



into Eq. (11)]; However, we have already seen evidence of anharmonic elastic effects in Ni- 
Au. Thus, we have performed LDA calculations of q(a±,k) directly from its definition in 
Eq. (Q), rather than using the harmonic approximation in Eq. (||). In Fig. CD, we show the 
results of the LAPW calculations of q Kl (a±,k) and q Au (a±,k) for six principle directions: 
(100), (111), (110), (201), (311), and (221). It is clear that the calculated values of q are 
not independent of a±, but rather show a marked and non-trivial dependence on the per- 
pendicular lattice constant. Thus, the lattice mismatch in Ni-Au appears to be too large 
for a harmonic continuum model of elasticity to be accurate. In particular, the value of 
g Nl (a_i_, 100) is quite low upon expansion, indicating that Ni is elastically extremely soft in 
this direction. Au, on the other hand, becomes softest in the (201) direction for significant 
compression. In a separate publication, [33] we will demonstrate that the anharmonic effects 



can be cast analytically in terms of the harmonic expressions of Eqs. (§)-@ by extending 
the expansion of ^(k): 

7(a ± , k) = a L {a ± )K L (k) (13) 

L 

to include angular momenta L=6,8, and 10 with the coefficients cll((i±) obtained from LDA 
calculations rather than the L=0,4 expression of Eq. @ used before. p3| 



The results for q Nl (a±, k) and q An (a±, k) are used to numerically minimize Eq. (|TT|) and 
hence to find AE^ s (k, x). The results for the CS energies are shown in Fig. ^| Here, also, 
the anharmonic effects are seen quite strongly as AE^? s (k, x) for some directions cross with 
other directions and asymmetries of the various directions are not all the same (effects which 
could not occur in the harmonic model). The most prominent feature of AE^? s (k, x) is that 
(100) is the softest elastic direction, which stems from the elastic softness of Ni along this 
direction. Ni being soft and Au being relatively hard along (100) leads to Ni(Au) being 
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highly distorted (nearly undistorted) for long-period (100) Ni-Au superlattices, and also 
leads to the marked asymmetry in Ai?Qg(100, x) towards the Ni-rich compositions. Similar 
arguments can be applied to explain the opposite asymmetry of the (201) strain. 

For E re { to be useful in the k-space CE, one must be able to know this energy for all 
directions, not merely the ones for which it was calculated. To obtain such a useful form, 
we fit the constituent strain results of Fig. ||] to a series of Kubic harmonics (0-10th order) 
consistent with cubic symmetry (L = 0,4,6,8,10). This procedure provides not only a good 
fit of the calculated strain data, but also an analytic form to obtain the values of AE^(k, x) 
for all directions. 

D. Stability of the cluster expansion 

Using the calculated formation energies {AH a } (Table PTI| ) and the anharmonic CS 
strain energy (Fig. |2[), we then fit the coefficients J(k) and {<//} of the k-space CE using 
Eq. (Q). We used all 33 calculated structures in the fit of the expansion, which included 20 
pair, 5 triplet, and 3 quadruplet interactions. The standard deviation of the fitted energies 
relative to their LAPW values is 5.3 meV/atom, which is the same order of magnitude as 
the numerical uncertainties in LAPW. The results for pair and multibody interactions are 
shown in Fig. |3|. 

In order for the expansion to have a useful predictive capability, tests must be performed 
to assess the stability of the fit: 

Changing the number of interactions: We performed tests of the stability of the fit with 
respect to the number of pair interactions, N paiis = (1 — 50). Figure [| shows the standard 
deviation of the fit as a function of the number of pairs interactions included. It is clear 
that the fit is well converged for N pa i TS = 20. We also tested the stability of the fit with 
respect to inclusion of more multibody interactions than are shown in Fig. ^: Including 
three additional triplet figures in the fit resulted in no change of the standard deviation 
of the fit, the added interactions had values < 2 meV/atom, and the original interactions 
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were changed by less than 2 meV/atom. Thus, the fit is stable with respect to the figures 
included (both pair and multibody). 

Changing the number of structures: We also performed tests of the predictive ability 
of the fit by removing some structures from the fit. First, we removed three structures 
which were originally fit quite well: Z2, (32, and Ll 2 (MAU3). Removing these structures 
from the input set resulted in their energies changing by <1 meV/atom. However, a much 
more critical test of the fit is to remove the structures which are fit most poorly: SQS14 a 
and SQS14b. Removing these structures from the fit changes their energies by only ~2- 
3 meV/atom. Thus, we are confident that the present k-space CE fit is both stable and 
predictive. 

IV. RESULTS OF CURRENT CALCULATIONS 

A. Mixing enthalpy: How good are previous calculations? 

Using the k-space cluster expansion in combination with a mixed real/reciprocal space 
Monte Carlo code (canonical), one can obtain thermodynamic properties of Ni-Au alloys. 
Figure |5] shows the mixing enthalpy as a function of temperature, Aif m ; x (T). Monte Carlo 
calculations were performed for a 16 3 =4096 atom cell, with 100 Monte Carlo steps per site 
for averages. The simulation was started at an extremely high temperature, and slowly 
cooled down using a simulated annealing algorithm. Also shown in Fig. [| is the value 
of the mixing energy of the completely random alloy. The difference between the Monte 
Carlo calculated Aif mix (T) and the random alloy energy is precisely the energetic effect of 
short-range order. We have fit the values of Aif m ; x (T) to linear and quadratic functions of 
(3 = 1/ksT to extrapolate the values down in temperature below the point at which coherent 
phase separation occurs in the simulation. (Both fits gave virtually identical results, so the 
linear fit is used here and below.) This allows us to ascertain the value of the mixing 
enthalpy at 1100 K, near the temperature where this quantity has been experimentally 
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measured. These results are tabulated in Table [IV], which shows both the effects of atomic 
relaxation (~ 100 meV/atom) and SRO (~ 25 meV/atom) on the mixing enthalpy, and 
compares the value of atomically relaxed and short-range ordered mixing energy with those 
values from experiment. One can see that by taking into account both relaxation and SRO, 
LDA produces a value for the mixing energy which is only different from experiment by 15- 
20 meV/atom. Thus, we conclude from this comparison that high quality LDA calculations 
provide accurate energetics for the Ni-Au system. 

The preceding discussion leads to a number of conclusions regarding previous calculations 
of AH mix : 

(i) Since relaxation reduces AH mix by ~100 meV/atom, the unrelaxed AH mix values ("d" 
in Table |) have to be reduced by this amount to appropriately compare with experiment. 

(ii) Since SRO reduces AH mix by ~25 meV/atom, the results of previous calculations 
that omitted SRO (all except "i" in Table |) have to be adjusted accordingly. 

(hi) In light of the fact that the empirical potential-based and ASA-based methods 
(LMTO and ASW) were shown to be inaccurate with respect to full-potential LDA methods 
for unrelaxed, ordered compounds (Table 0), the results of relaxed, mixing energies of random 
alloys appear to be questionable using these schemes. 

B. Configurational or non-configurational entropy? 

From the fit of the Monte Carlo data in Fig. ^, one can find the configurational entropy 
of the Ni .5Auo.5 disordered phase by integrating the energy down from infinite temperature 
(where the configurational entropy is known): 

AS(T) = AS(T = oo) + E(T)/T - k B f E((3)d(3 (14) 

Jo 

The configurational entropy obtained from thermodynamic integration in this way is 

A5 conf . (Ni . 5 Au .5, T = HOOK) = 0.56fc B , (15) 
compared to the "ideal" (infinite temperature) value of 
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AS conf .(Ni . 5 Au .5,T -> oo) = 0.69k B , (16) 

This calculated value for the configurational entropy of mixing can be compared with the 
experimentally measured values of total entropy of mixing: Calorimetric measurements give 
AS(T = 1150/0 = Q while EMF measurements give AS(T = 1173iT) = 1.08Jfc B 

||1 1|| . Thus, we can obtain an estimate of the non- configurational entropy, and find it to be 
large: AS^on-conf . (T ~ HOOK) = 1.04 — 0.56 = 0A8k B . This non-configurational entropy 
is hence responsible for Tmg being so small experimentally, compared to all the theoret- 
ical results. In fact, if we use the calculated Aif m ; x = 93 meV/atom and the combined 
"experimental/calculated" AS^on-conf. = 0.48fcg in the following formula: 



'-MG 



(17) 



k B + 2AS' non 

—conf. 

we obtain Tmg ~ 1100 K and &bTmg/ Af/ mix = 1.02, much closer to the experimental 
values (Tmg ~ 1083 K and k B T MG / AH mix = 1.2) than using the above formula neglecting 
non-configurational entropy (Tmg ~ 2150 K and k B TMc/ AH mix = 2.0). 

From this consideration of non-configurational effects, one should conclude that the accu- 
racy of a calculation with configurational degrees of freedom only (as is done in most of the 
previous calculations |]3"5|), should be determined by looking at the energetics, not the tran- 
sition temperatures. Thus, previous calculations which give "good" transition temperatures 
do so precisely because they have "bad" energetics. 



C. Short-range order of Nii_ z Au x solid solutions 

Using the k-space CE and Monte Carlo, we may also compute the SRO of disordered 
Nii_ x Au a; alloys. We show the results of our SRO simulations for Ni .4Auo.6 in Fig. |6]. For 
the SRO Monte Carlo calculations, a cell of 24 3 =13824 atoms was used, with 100 Monte 
Carlo steps for equilibration, with averages taken over the subsequent 500 steps. Several 
calculations and measurements of the SRO exist in the literature: Wu and Cohen used 
diffuse x-ray scattering to deduce the atomic SRO of Ni .4Au .6 at T=1023 K. The measured 
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diffuse intensity due to SRO must be separated from all the other contributions which give 
rise to diffuse intensity, and for this purpose, Wu and Cohen used 25 real-space Fourier 
shells of SRO parameters, and found the rather surprising result that the peak intensity 
in reciprocal space due to SRO is of ordering-type and occurs at the point ksRo=(0.6,0,0), 
rather than k SRO =(0,0,0) which would be expected for a clustering alloy. Several authors 
have tried to account for this ordering nature of the SRO: Lu and Zunger || calculated the 
SRO (using 21 real-space shells) and found peaks at ~(0.8,0,0) whereas Asta and Foiles [17| 



used an embedded atom method and found the SRO (using 8 real-space shells) to peak at 
~(0.5,0,0). Our calculations for the SRO of M0.4AU0.6 are given in Fig. || We have calculated 
the SRO at T=2300 K, above the miscibility gap temperature for our alloy Hamiltonian. 
We find that, using 8, 25, and 100 shells, the SRO peaks at (0.65,0,0), (0.40,0,0), and 
(0.38,0,0) respectively, in good agreement with both the measurements of Wu and Cohen 
[ksRo=(0.6,0,0) for 25 shells] and also with previous calculations. 

Equation (|1^) shows that the alloy Hamiltonian used in the Monte Carlo simulations is 
composed of three parts: the pair interaction terms, the multibody interaction terms, and 
the constituent strain terms. It is interesting to see the effect of each of these portions of the 
alloy Hamiltonian on SRO. Thus, in addition to the "full" calculations, which contain pairs, 
multibodies, and constituent strain in the alloy Hamiltonian, we have also computed the 
SRO with (i) the CS energy only, and (ii) the CS energy plus the pair interactions. These 
results are shown in Fig. |7|. (Because the CS energy is non-analytic in reciprocal space about 
the origin, many Fourier coefficients are required to converge the SRO of CS alone, thus we 
show only results using 100 shells of parameters in Fig. |7|.) One can see that the SRO with 
CS only is dominated by almost constant streaks of intensity along the T — X line, and very 
little intensity elsewhere. This SRO pattern is understandable when one considers that the 
constituent strain at this composition (Fig. is much softer (much lower in energy) in the 
(100) direction than along any other direction. Thus, (100)-type fluctuations in the random 
alloy are be energetically favored, and because the constituent strain is dependent only on 
direction and not on the length of the wavevector, one should expect that all fluctuations 
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along the (100) will occur roughly equally, regardless of the length of the wavevector. This 
is precisely what we see in Fig. [T[ Contrasting this SRO using CS only with that calculated 
both CS energy and pair interactions (but not multibody interactions) shows that the pair 
interactions create a peak in intensity along the T — X line, but significantly closer to V 
than the peak intensity using the "full" alloy Hamiltonian. Thus, while the effect of pairs 
is to create a peak near the T point, the multibody interactions move this peak out from T 
towards the X-point. 



D. Standard inverse Monte Carlo would give unphysical interaction energies: a 

challenge 

The statistical problem we have solved here involves the calculation of the alloy SRO 
at high temperature for given alloy Hamiltonian ({Jij}, {<//}, and AEqs). However, a 
popular technique used to study phase stability in alloys involves the "inverse" problem of 
determining a set of pair-only interactions {Jij} from a measured or calculated SRO pattern, 
and the subsequent use of these pair interactions to determine thermodynamic properties 
other than the SRO. In fact, {Jij} are often used to determine AiJ mix or phase stability. 
As we have mentioned in the introduction and described more fully in Ref. |3q| , inverting 



the SRO always removes information on energy terms that are SRO-independent, e.g., the 
volume deformation energy G(x). This loss prevents, in principle, the interactions deduced 
from SRO from being applied to predict physical properties which depend on G(x), such as 
A// mix . For example, in the case of Ni-Au, the SRO is of ordering-type. Thus, we expect 
that inverting the SRO of Ni-Au (e.g., via inverse Monte Carlo) would produce interactions 
{Jij} which are of ordering-type, and using these interactions to predict the mixing enthalpy 
would result in the unphysical result AH mix < 0. 

One might suspect that by changing the temperature, one could obtain a shift of the 
SRO from ordering- to clustering-type, and thus, the inverse technique would then produce 
interactions which would correctly give AH mix > 0. However, we have computed the SRO 
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for several temperatures, and find no evidence of a shift in SRO to clustering-type. 

A test of our expectations by any of the practitioners of inverse Monte Carlo would 
certainly be welcomed. To that end, our SRO calculations are available for use as input 
to inverse Monte Carlo to extract interactions. These SRO calculations are available for 
a variety of compositions and temperatures, each with a large number of real-space SRO 
parameters. It would be of great interest to see whether the interactions extracted from 
inverting the SRO of Ni-Au would produce the correct sign of AH m[x . 
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FIGURES 

FIG. 1. LAPW calculations of q W (a ± , k) of Eq. (|) for Ni-Au. Shown are (a) q m and (b) q 
for six principle directions. 



FIG. 2. LAPW calculations of AEcs(k,x) for Ni-Au for six principle directions. 

FIG. 3. (a) Pair and (b) multibody interaction energies for Ni-Au. The multibody figures 
are defined by the following lattice sites, in units of a=2 (the origin is contained in all fig- 
ures): J 3 - (110),(101), K 3 - (110),(200), N 3 - (200), (002), P 3 - (110),(103), Q 3 - (110), (220), 
J 4 - (110), (101), (011), K 4 - (110), (101), (200), and L 4 - (110), (101), (211). 

FIG. 4. Cluster expansion fitting error in Ni-Au versus the number of pair interactions included 
in the fit. 

FIG. 5. AH(T) computed for Nio.sAuo.s from a combination of the k-space cluster expansion 
and Monte Carlo simulations. 

FIG. 6. Monte Carlo-calculated short-range order of Nio.4Auo.6 in the (hkO) plane using (a) 8, 
(b) 25, and (c) 100 shells of Warren Cowley SRO parameters. Peak intensity is red shaded contour 
while the lowest contours are shaded blue. Contours are separated by 0.1 Laue unit in each plot. 

FIG. 7. Monte Carlo-calculated short-range order of M0.4AU0.6 using (a) constituent strain 
terms only, (b) constituent strain and pair terms, and (c) constituent strain, pair, and multibody 
terms in the alloy Hamiltonian. Peak intensity is red shaded contour while the lowest contours are 
shaded blue. Contours are separated by 0.1 Laue unit in each plot. 
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TABLES 

TABLE I. Summary of energy calculations performed for Nii_ x Au x alloys. Shown are the 
methods used to compute T=0 energetics, as well as the type of cluster expansion (CE) and 
statistics used. Also given is the mixing energy of the T — > oo random alloy near x=l/2, and the 
calculated miscibility gap temperature, if available. FLAPW = full-potential linearized augmented 
plane wave method, FLMTO = full-potential linear muffin-tin orbital method, ASW = augmented 
spherical wave method, LMTO-ASA = linear muffin-tin orbital method in the atomic sphere ap- 
proximation, EAM = embedded atom method, MC = Monte Carlo, CVM = cluster variation 
method, MF = mean-field, SOE = second-order expansion. 



Authors 


T=0 Energy 


Method 

Cluster 
Expansion 
Technique 


Statistics 


Results 


Wolverton and Zunger a 


FLAPW 


k-space CE 


MC 


+118 




Lu and Zunger 


FLAPW 


t-G 


MC 


+127 




Dcutsch and Pasturel c 


FLMTO 


€ - G 


none 


+136 




Takizawa, Terakura, and Mohri d 


ASW 


CW 


CVM 


+170 




Amador and Bozzolo e 


LMTO-ASA 


cw 


CVM 


+150 




Colinet et al? 


LMTO-ASA 


t-G 


CVM 


+67 


1200-1400 


Morgan and de Fontaine s 


LMTO-ASA + 


e-G 


CVM 


+98 


2330 




"Elastic Springs" 










Eymery et al. 


Empir. Potential 


Simulation 


none 


+60 




Tetot and Finel 1 


Empir. Potential 


Simulation 


MC 


+48 m 


950 


Dcutsch and Pasturel c 


Empir. Potential 


Simulation 


none 


+83 




Asta and Foiles J 


EAM 


SOE 


MC/MF 


+78 


2460 


Expt. (Calorimctry) T=1150 K k 
Expt. (EMF) T=1173K 1 
Expt. (Phase Diagram) 




+76 
+77 


1083 



a Present results. 
b Ref. | 



=Ref. 112 



d Ref. p| 



e Ref. |3| 
f Ref. 
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SRef. p| 
h Ref. |9| 
iRef. || 
JRef. |7| 
k Ref. § 
! Ref. Q 
m at T=1150 K 



TABLE II. Comparison of formation enthalpies AHf(a) for Ni-Au ordered compounds. 



Nomenclature for the ordered structures is the same as that used in Ref. \22\. All energies in 
meV/atom. Numbers in parentheses indicate unrelaxed energies. 





Fully Relaxed 


Pctrticilly Relaxed 








Empirical 




Structure 


FLAPW a 


FLAPW b 


FLMTO c 


ASW d 


LMTO 6 


Potential c 


EAM J 


NiAu (Ll ) 


+76.1 (+98.1) 


+76.8 


+79.4 (+96.4) 


(+59) 


(+116.6) 


+57.9(+73.9) 


+21.4(4-91.1) 


Ni 2 Au 2 (Z2) 


+70.2 (+286.7) 


+124.3 


+123.1 (+300.1) 




(+213.4) 


+62.3(+127.7) 


-130.3(4208.6) 


NiAu (ili) 


+166.8 (+192.3) 


+167.6 


+175.4 




(+177.9) 




+72.9(+159.7) 


NiAu ("40") 


+84.8 (+93.5) 


+83.8 


+89.9 




(+114.3) 




-1.9(496.4) 


Ni 3 Au (Ll 2 ) 


+77.5 


+75.5 


+80.7 


+42 


+92.4 


+58.1 


477.1 


Ni 3 Au (DO22) 


+75.0 (+75.0) 




+81.5 




(+95.3) 






NiAu 3 (Ll 2 ) 


+78.9 


+78.2 


+78.0 


+52 


+89.4 


+54.1 


486.1 


NiAu 3 (DO22) 


+68.6 (+68.7) 




+68.0 




(+76.4) 







a Present results. Complete atomic relaxation via quantum mechanical forces and total-energy 
minimization. 

b Ref. P]. Partial atomic relaxation via continuum elasticity, using Eqs. 



c Ref. g2| 
d Ref. p| 



c Ref. [15]. LMTO- AS A with sphere radii chosen to minimize charge transfer. 
jRef. |7| 
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TABLE III. Listing of the LAPW calculated unrelaxed and relaxed AH (a) [in meV/atom] for 
Nii-^Aoj;. Many of the structures calculated here can be characterized as a (Ni) p (Au) g superlattice 



of orientation k. We use the nomenclature of Ref. |2^| for structure names. 



Orientation 


fooil 


foul 


10121 


fllll 

L J 


[1131 


triTTYl 1 1 1 a 
1 VJ1 1 1 1 LI 1 cX 












AB 


Ll 


Ll 


Ll 


Lli 


Lli 


Unrelaxed 


+98.1 


+98.1 


+98.1 


+192.3 


+192.3 


R claxcd 


+76.1 


+76.1 


+76.1 


+166.8 


+166.8 


CE( Relaxed) 


+74.8 


+74.8 


+74.8 


+167.1 


+167.1 


A 2 B 


31 


t1 
1 x 




al 




Unrelaxed 


+207.8 


+123.3 




+288.5 




H claxcd 


+105.7 


+98.9 




+202.2 




CE( Relaxed) 


+105.9 


+102.4 




+208.4 




AB2 








a 2 




Unrelaxed 


+151.7 


+126.3 




+200.9 




claxcd 


+38.3 


+102.6 




+100.9 




CF^RelaxpH") 

V > 1 J I 1_ tCldACU j 


+37.8 


+98.8 




+94.5 




A 3 B 


Zl 


Yl 


DO22 


VI 


Wl 


Unrelaxed 


+221.7 


+148.5 


+75.0 


+290.8 




Relaxed 


+89.9 


+99.2 


+75.0 


+193.7 


+125.7 


CE( Relaxed) 


+94.3 


+91.3 


+69.1 


+189.6 


+120.8 


AB3 


Z3 


Y3 


DO22 


V3 


W3 


Unrelaxed 


+142.0 


+104.1 


+68.7 


+172.8 




Relaxed 


+32.4 


+78.7 


+68.6 


+83.0 


+88.4 


CE( Relaxed) 


+28.2 


+77.7 


+67.6 


+79.1 


+83.2 


A2-B2 


Z2 


Y2 


"40" 


V2 


W2 


Unrelaxed 


+286.7 


+192.3 


+93.5 


+335.8 


+144.2 


Relaxed 


+70.2 


+96.6 


+84.8 


+162.4 


+93.6 


CE(Relaxed) 


+69.9 


+101.1 


+88.3 


+166.7 


+99.3 














Unrelaxed 


+576.2 


+576.2 


+576.2 


+576.2 


+576.2 


Relaxed 


+30.8 


+117.7 


+84.8 


+173.8 


+119.8 


CE(Relaxed) 


+30.8 


+116.1 


+86.8 


+172.5 


+117.9 




Other 


Structures 










LI2 (A 3 B) 


Ll 2 (AB 3 ) 


D7 (A7B) 


D7 b (A 7 B) 




Unrelaxed 


+77.5 


+78.9 


+82.9 


56.8 




Relaxed 


+77.5 


+78.9 


+82.9 


56.8 




CE(Relaxed) 


+80.7 


+78.6 


+98.5 


57.6 






SQS14 a (A 6 B 2 ) 


SQS14i, (A 2 B 6 ) 


Z6 (A3B3 - 100) 


Z5 (A 2 B 3 - 100) 




Unrelaxed 


+183.2 


+118.2 


+355.5 


+273.3 
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Relaxed +96.8 +59.8 +63.2 +57.1 

CE(Relaxed) +81.5 +75.0 +62.5 +57.9 
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TABLE IV. Ai7 m ix for Nio.5Auo.5- All energies in meV/atom. SQS-4 refers to a 4-atom special 
quasi-random structure (Y2). This table shows the effects of relaxation (first line minus second 
line) and short-range order (third line minus fourth line) on the mixing energy. 



AH ■ 


SQS-4 Unrelaxed(T = 00) 


+ 192 


SQS-4 Relaxed(T = 00) 


+97 


CE Relaxed (T = 00) 


+ 118 


CE Relaxed(T = HOOK) 


+93 


Expt. (Calorimetry) T=1150 K 


+76 


Expt. (EMF) T=1175 K 


+77 
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6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 
a, (a.u.) 




6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 
a, (a.u.) 
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Ni-Au Pair Interactions 
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(b) 




Ni-Au Multibody Interactions 




3-body 






4-body 
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